function [ out ] = osl_tf_transform( S, dat )
%
% [ dattf ] = osl_tf_transform( S , dat )
%
% Either downsamples the time domain data ('S.tf_method = timedomain_ds')
% or does a time-frequency transformation.
%
% If called with only S, and no dat, this function will only produce output
% containers and not do the transform.  If passed with a second argument
% dat, it will do the transformation on the data in dat.  dat should have
% dimensions [nTrials,nTimepoints]
%
% ARGUMENTS (in struct 'S')
% - S.tf_method                 : e.g. 'none', 'morlet', 'hilbert', 'multitaper'
% - S.tf_freq_range             : the frequency window over which the tf is performed, e.g. [4 35]
% - S.tf_num_freqs                 : how many frequency bands to compute?
% - S.raw_times                  : the times in source_recon_results.times
% - S.tf_logtransform           : log transform the tf, trial by trial?
% - S.ds_factor                 : used for time domain, hilbert, morlet
% - S.tf_morlet_factor             : only used for morlet
% - S.tf_multitaper_twin                   : used for multitaper
% - S.tf_multitaper_taper                  : if multitaper, hanning or dpss
% - S.tf_multitaper_freqsmooth             : used for multitaper
% - S.tf_calc_amplitude
%
% IF S.ft_method = 'morlet', optional argument:
% - S.tf_morlet_basis (create using 'spm_eeg_morlet', for e.g.)
% >NB if this is NOT set, it will be reset each time the function is called
%
% OUTPUTS (in struct 'out')
% - out.tf_times
% - out.tf_freqs
% - out.tf_freq_res
%
% IF S.doTransform = 1 then additionally
% - out.dattf
%
% IF doing the HILBERT, then additionally
% - out.datbp
%
% IF doing the MORLET, then additionally
% - out.tf_morlet_basis

if nargin < 2
    S.doTransform = 0;
    dat = [];
else
    S.doTransform = 1;
end % if nargin < 2

if(~isfield(S,'tf_calc_amplitude')),
    S.tf_calc_amplitude=1;
end;

% Do downsampling / TF transform
switch S.tf_method
    
    case 'none'
        
        % always do
        out.tf_freqs = -1; % time domain
        if ~isempty(S.ds_factor)
            [dsn dsd]    = rat(S.ds_factor);
            tmp          = resample(S.raw_times,dsn,dsd);
            out.tf_times = linspace(S.raw_times(1),S.raw_times(end),numel(tmp));
            clear tmp
        else
            out.tf_times  = S.raw_times;
        end
        
        % Continue and do the transform if S.doTransform == 1
        if S.doTransform
            
            if ~isempty(S.ds_factor)
                [dsn dsd]    = rat(S.ds_factor);
                dat = dat';
                out.dattf    = resample(dat,dsn,dsd);
                out.dattf = out.dattf';
            else
                out.dattf = dat;
            end % if ~isempty(S.ds_factor)
        else
            
            % the first time the function is called, print some commentary
            if isempty(S.ds_factor)
                tdHz  = 1 ./ (S.raw_times(2) - S.raw_times(1));
                glmHz = 1 ./ (out.tf_times(2) - out.tf_times(1));
                disp(['TIME DOMAIN, NOT DOWNSAMPLING; raw data at ' num2str(tdHz) 'Hz, GLM at ' num2str(glmHz) 'Hz.'])
            else
                tdHz  = 1 ./ (S.raw_times(2) - S.raw_times(1));
                glmHz = 1 ./ (out.tf_times(2) - out.tf_times(1));
                disp(['TIME DOMAIN, DOWNSAMPLE FACTOR = ' num2str(S.ds_factor) '; raw data at ' num2str(tdHz) 'Hz, GLM at ' num2str(glmHz) 'Hz.']);
            end % if isempty(S.time_downsample_factor)
            
        end % if S.doTransform
        
    case 'hilbert'
        
        % always do
        if S.tf_num_freqs > 1
            out.tf_freq_res = S.tf_hilbert_freq_res;
            out.tf_freqs   = linspace(S.tf_freq_range(1)+out.tf_freq_res/2, S.tf_freq_range(2)-out.tf_freq_res/2, S.tf_num_freqs);
        else
            out.tf_freq_res = S.tf_hilbert_freq_res;
            out.tf_freqs   = (S.tf_freq_range(2) + S.tf_freq_range(1)) ./ 2;
        end % if S.tf_num_freqs > 1
        
        % downsample the time vector appropriately
        if ~isempty(S.ds_factor)
            [dsn dsd]    = rat(S.ds_factor);
            tmp          = resample(S.raw_times,dsn,dsd);
            out.tf_times = linspace(S.raw_times(1),S.raw_times(end),numel(tmp));
            clear tmp
        else
            out.tf_times  = S.raw_times;
        end
        
        % Continue and do the transform if S.doTransform == 1
        if S.doTransform
            
            tempoutpow = zeros(size(dat,1),length(S.raw_times),length(out.tf_freqs));
            tempoutbp = zeros(size(dat,1),length(S.raw_times),length(out.tf_freqs));
            fsample=1/(S.raw_times(2)-S.raw_times(1));
            out.bp_cutoffs=zeros(length(out.tf_freqs),2);
            
            for iTrial = 1:size(dat,1); % indexes trials
                
                dattrl=dat(iTrial,:);
                dattrl=dattrl(:);
                
                for f = 1 : length(out.tf_freqs)
                    
                    % bandpass data
                    bp_cutoff=max(1e-3,[out.tf_freqs(f)-out.tf_freq_res, out.tf_freqs(f)+out.tf_freq_res]);
                    tempdatbp=bandpass(dattrl,bp_cutoff,fsample,0);
                    %tempdatbp = ft_preproc_bandpassfilter(dattrl, fsample, bp_cutoff, 3);
                    out.bp_cutoffs(f,:)=bp_cutoff;
                    
                    signal_h = hilbert(tempdatbp);
                    
                    % square the hilbert transform and take the
                    % +ve square root to get the estimated envelope
                    if(S.tf_calc_amplitude)
                        tempdatpow = (sqrt(signal_h.*conj(signal_h)));
                    else
                        tempdatpow=signal_h;
                    end;
                    
                    tempoutpow(iTrial,:,f) = tempdatpow;
                    tempoutbp(iTrial,:,f)  = tempdatbp;
                    
                end % f = 1 : length(out.tf_freqs)
            end % for iTrial = 1:size(dat,1);
            
            % if requested, downsample
            if ~isempty(S.ds_factor)
                [dsn dsd]    = rat(S.ds_factor);
                
                nfreqs = size(tempoutpow,3);
                
                tempoutpow   = permute(tempoutpow,[2,1,3]);
                tempoutbp    = permute(tempoutbp,[2,1,3]);
                
                tempoutpow   = reshape(tempoutpow,[size(tempoutpow,1),(size(tempoutpow,2)*nfreqs)]);
                tempoutbp    = reshape(tempoutbp, [size(tempoutbp,1), (size(tempoutbp, 2)*nfreqs)]);

                tempoutpow   = resample(tempoutpow,dsn,dsd);
                tempoutbp    = resample(tempoutbp,dsn,dsd);
                
                tempoutpow   = reshape(tempoutpow,[size(tempoutpow,1),(size(tempoutpow,2)./nfreqs),nfreqs]);
                tempoutbp    = reshape(tempoutbp,[size(tempoutbp,1), (size(tempoutbp,2)./nfreqs),nfreqs]);
              
                out.dattf = permute(tempoutpow,[2,1,3]);
                clear tempoutpow
                out.datbp = permute(tempoutbp,[2,1,3]);
                clear tempoutbp

            else
                out.dattf   = tempoutpow;
                clear tempoutpow
                out.datbp   = tempoutbp;
                clear tempoutbp
            end
            
            if S.tf_logtransform
                out.dattf = log(out.dattf);
            end % S.tf_logtransform
            
        else
            
            % the first time the function is called, print some commentary
            disp('USING HILBERT TF');
            
            disp(['Freqs: ' num2str(out.tf_freqs)]);
            disp(['Freq resolution: ' num2str(out.tf_freq_res)]);
            
        end % if S.doTransform
        
        
    case 'morlet'
        
        % always do
        if S.tf_num_freqs == 1
            out.tf_freqs   = (S.tf_freq_range(1) + S.tf_freq_range(2))/2;
            out.tf_freq_res =  S.tf_freq_range(1) - S.tf_freq_range(2);
        else
            out.tf_freqs   = linspace(S.tf_freq_range(1),S.tf_freq_range(2),S.tf_num_freqs);
            out.tf_freq_res = out.tf_freqs(2) - out.tf_freqs(1);
        end % if S.tf_num_freqs == 1
        
        % downsample the time vector appropriately
        if ~isempty(S.ds_factor)
            [dsn dsd]    = rat(S.ds_factor);
            tmp          = resample(S.raw_times,dsn,dsd);
            out.tf_times = linspace(S.raw_times(1),S.raw_times(end),numel(tmp));
            clear tmp
        else
            out.tf_times  = S.raw_times;
        end % if ~isempty(S.ds_factor)
        
        if ~isfield(S,'tf_morlet_basis')
            disp('Creating Morlet basis set.  If you are seeing message many times, you may wish to pass a morelt basis set to ''osl_tf_transform''');
            fres=1./(S.raw_times(2)-S.raw_times(1));
            out.tf_morlet_basis = spm_eeg_morlet(S.tf_morlet_factor, 1000/fres, out.tf_freqs); % ?? GW why is sample time coded this way?
            S.tf_morlet_basis   = out.tf_morlet_basis;
        end % if ~isfield(S,'tf_morlet_basis')
        
        % Continue and do the transform if S.doTransform == 1
        if S.doTransform
            
            tempdattf = zeros(size(dat,1),length(S.raw_times), length(out.tf_freqs));
            for iTrial = 1:size(dat,1); % indexes trials
                
                dattrl = dat(iTrial,:);
                dattrl = dattrl(:);
                
                for f = 1 : length(out.tf_freqs)
                    tmp = conv(dattrl, S.tf_morlet_basis{f},'same');
                    
                    % power
                    if(S.tf_calc_amplitude)
                        tmp = sqrt(tmp.*conj(tmp))'; % tf_num_freqs x num_timepoint
                    end;
                    
                    tempdattf(iTrial,:,f) = tmp;
                end % for f = 1 : length(out.tf_freqs)
            end % for iTrial = 1:size(dat,1);
            
            % downsample
            if ~isempty(S.ds_factor)
                [dsn dsd]    = rat(S.ds_factor);
                
                nfreqs = size(tempdattf,3);
                
                tempdattf   = permute(tempdattf,[2,1,3]);
                tempdattf   = reshape(tempdattf,[size(tempdattf,1),(size(tempdattf,2)*nfreqs)]);
                tempdattf   = resample(tempdattf,dsn,dsd);
                tempdattf   = reshape(tempdattf,[size(tempdattf,1),(size(tempdattf,2)./nfreqs),nfreqs]);
              
                out.dattf = permute(tempdattf,[2,1,3]);
                clear tempdattf;
            else
                out.dattf = tempdattf;
            end
            
            if S.tf_logtransform
                out.dattf(:,f) = log(out.dattf(:,f));
            end % if S.tf_logtransform
        else
            
            % the first time the function is called, print some commentary
            disp('USING MORLET TF');
            
            disp(['Freqs: ' num2str(out.tf_freqs)]);
            
        end % if S.doTransform
        
        
        
    case 'multitaper'
        
        if S.tf_multitaper_twin < (1/S.tf_freq_range(1))
            error('The maximum time window specified for the analysis is less than one cycle of the lowest frequency!')
        end
        
        % always do
        
        switch S.tf_multitaper_taper
            
            case 'hanning'
                
                % always do
                
                % DETERMINE THE FREQUENCY AXIS
                
                % We need to set the frequency bands such that we have a full number
                % of cycles within the time window.
                out.tf_freqs   = linspace(S.tf_freq_range(1),S.tf_freq_range(2),S.tf_num_freqs);
                cycles = S.tf_multitaper_twin.*out.tf_freqs;
                fullCycles = floor(cycles);
                if ~isempty(S.tf_multitaper_ncycles)
                    fullCycles = min(fullCycles,S.tf_multitaper_ncycles);
                end
                
                % DETERMINE THE TIME AXIS
                
                % determine the true time range we can use
                first_viable_time = S.raw_times(1)   + S.tf_multitaper_twin/2;
                last_viable_time  = S.raw_times(end) - S.tf_multitaper_twin/2;
                
                rawTstep        = S.raw_times(2) - S.raw_times(1);
                if isempty(S.ds_factor)
                    S.ds_factor = 0.1;
                end % if isempty(S.ds_factor)
                dsTstep         = rawTstep./S.ds_factor;
                out.tf_times    = first_viable_time:dsTstep:last_viable_time;
                
                if S.doTransform % ... continue and do the transform
                    
                    fsample=1/(S.raw_times(2)-S.raw_times(1));
                    
                    % parse data matrix into fieldtrip format: cell array
                    % of trials; as this is voxelwise, each array will be
                    % {1xn}
                    ftraw = [];
                    ftraw.trial = cell(size(dat,1),1);
                    for iTrial = 1:size(dat,1);
                        ftraw.trial{iTrial} = dat(iTrial,:);
                    end % for iTrial = 1:size(dat,1);
                    ftraw.time = {S.raw_times};
                    ftraw.time = repmat(ftraw.time,[size(dat,1),1]);
                    ftraw.fsample = fsample;
                    ftraw.label{1, 1} = ['Ch'];
                    
                    cfg = [];
                    cfg.output      ='pow'; % NB could be 'fourier' - may need care though...
                    cfg.taper       ='hanning';
                    cfg.method      ='mtmconvol';
                    cfg.foi         = out.tf_freqs;
                    cfg.toi         = out.tf_times;
                    cfg.t_ftimwin   = fullCycles./cfg.foi;
                    cfg.keeptrials = 'yes';
                    
                    
                    freq = ft_freqanalysis(cfg, ftraw);
                    pow = freq.powspctrm;
                    
                    % in the first level function, dat is [trials,times,freqs]
                    out.dattf = permute(pow,[1 4 3 2]); % this permute operation gets rid of the channel dimension, as well as rearranging the order of the other dimensions
                    
                    
                    %                     % deal with the special cases where nTimes or nFreqs =
                    %                     % 1;
                    %                     if (numel(cfg.toi) == 1) && (numel(cfg.foi) == 1) % a single tf value!
                    %
                    %                         % out.dattf =
                    %
                    %                     elseif numel(cfg.toi) == 1; % only one time
                    %
                    %
                    %
                    %                     elseif numel(cfg.foi) == 1; % only one frequency
                    %
                    %                         % in the first level function, dat is [trials,times,freqs]
                    %                         out.dattf = permute(pow,[1 4 3 2]); % this permute operation gets rid of the channel dimension, as well as rearranging the order of the other dimensions
                    %
                    %
                    %                     else % multiple times and freqs
                    %
                    %                         % in the first level function, dat is [trials,times,freqs]
                    %                         out.dattf = permute(pow,[1 4 3 2]); % this permute operation gets rid of the channel dimension, as well as rearranging the order of the other dimensions
                    %
                    %                     end % if (numel(cfg.toi) == 1) && (numel(cfg.foi) == 1)
                    
                    if S.tf_logtransform % take the log
                        out.dattf = log(out.dattf);
                    end
                    
                else
                    
                    disp(['Hanning taper TF.']);
                    disp(['The first time for which the TF will be calculated will be ' num2str(first_viable_time) 's.']);
                    disp(['The last time for which the TF will be calculated will be ' num2str(out.tf_times(end)) 's.']);
                    
                end % if S.doTransform
                
                
            case 'dpss'
                
                % always do
                
                % DETERMINE THE FREQUENCY AXIS
                
                % We need to set the frequency bands such that we have a full number
                % of cycles within the time window.
                out.tf_freqs   = linspace(S.tf_freq_range(1),S.tf_freq_range(2),S.tf_num_freqs);
                cycles = S.tf_multitaper_twin.*out.tf_freqs;
                fullCycles = floor(cycles);
                if ~isempty(S.tf_multitaper_ncycles)
                    fullCycles = min(fullCycles,S.tf_multitaper_ncycles);
                end
                
                % DETERMINE THE TIME AXIS
                
                % determine the true time range we can use
                first_viable_time = S.raw_times(1)   + S.tf_multitaper_twin/2;
                last_viable_time  = S.raw_times(end) - S.tf_multitaper_twin/2;
                
                rawTstep        = S.raw_times(2) - S.raw_times(1);
                if isempty(S.ds_factor)
                    S.ds_factor = 0.1;
                end % if isempty(S.ds_factor)
                dsTstep         = rawTstep./S.ds_factor;
                out.tf_times    = first_viable_time:dsTstep:last_viable_time;
                
                if S.doTransform % ... continue and do the transform
                    
                    fsample=1/(S.raw_times(2)-S.raw_times(1));
                    
                    % parse data matrix into fieldtrip format: cell array
                    % of trials; as this is voxelwise, each array will be
                    % {1xn}
                    ftraw = [];
                    ftraw.trial = cell(size(dat,1),1);
                    for iTrial = 1:size(dat,1);
                        ftraw.trial{iTrial} = dat(iTrial,:);
                    end % for iTrial = 1:size(dat,1);
                    ftraw.time = {S.raw_times};
                    ftraw.time = repmat(ftraw.time,[size(dat,1),1]);
                    ftraw.fsample = fsample;
                    ftraw.label{1, 1} = ['Ch'];
                    
                    cfg = [];
                    cfg.output      ='pow'; % NB could be 'fourier' - may need care though...
                    cfg.taper       ='dpss';
                    cfg.method      ='mtmconvol';
                    cfg.foi         = out.tf_freqs;
                    cfg.toi         = out.tf_times;
                    cfg.t_ftimwin   = fullCycles./cfg.foi;
                    cfg.tapsmofrq   = S.tf_multitaper_freqsmooth.*cfg.foi;
                    cfg.keeptrials = 'yes';
                    
                    freq = ft_freqanalysis(cfg, ftraw);
                    pow = freq.powspctrm;
                    
                    % in the first level function, dat is [trials,times,freqs]
                    out.dattf = permute(pow,[1 4 3 2]); % this permute operation gets rid of the channel dimension, as well as rearranging the order of the other dimensions
                    
                    if S.tf_logtransform % take the log
                        out.dattf = log(out.dattf);
                    end
                    
                else
                    
                    disp(['Slepian taper TF.']);
                    disp(['The first time for which the TF will be calculated will be ' num2str(first_viable_time) 's.']);
                    disp(['The last time for which the TF will be calculated will be ' num2str(out.tf_times(end)) 's.']);
                    
                    
                end % if S.doTransform
                
        end % switch cfg.taper
        
    otherwise
        
        error('Invalid TF method');
        
end % switch S.tf_method